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. We present an analytical approach to deal with nonlinear delay differential 

equations close to instabilities of time periodic reference states. To this end 
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we start with approximately determining such reference states by extending 
the Poincare Lindstedt and the Shohat expansions which were originally de- 
veloped for ordinary differential equations. Then we systematically elaborate 
a linear stability analysis around a time periodic reference state. This allows 
to approximately calculate the Floquet eigenvalues and their corresponding 
eigensolutions by using matrix valued continued fractions. 
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I. INTRODUCTION 



Over the last two decades considerable new interest in the theory of delay differential 
equations led to various remarkable results [|r||. The reason is that the solution space for 
delay differential equations has to be considered as infinite dimensional although only a finite 
number of dynamical variables is involved 0. As a consequence, nonlinear delay differential 
equations reveal a broad class of instabilities leading from oscillatory to chaotic behavior. 
Apart from the period doubling route to chaos also quasi periodic states, intermittency and 
^ locking behavior have been observed in detailed numerical studies @. In particular in the 

chaotic domain it has been suggested that the envelope to the Kaplan Yorke dimension of 
a delay induced chaotic attractor is proportional to the time delay . This fact offers 

the possibility to generate high dimensional chaotic attractors by simply increasing the time 
delay. 

Delay differential equations have been successfully applied to model numerous nonlinear 
systems where dynamical instabilities are induced by the finite propagation time of signals 
in feedback loops. For instance, experiments on optical devices, acusto optic and electro 



optic bistable devices |I]]-4U| have confirmed both the theoretical and the numerical predic- 
tions. But delay induced instabilities play also an important role in other disciplines such as 
population dynamics ||, radio engineering sciences [14H, economy [15 and biology [16|. In 



addition it has been noted in medical sciences that there exists a remarkable variety of clin- 
ically relevant dynamical phenomena under physiological and pathological conditions. For 
example oscillations or chaotic behavior can spontaneously occur or disappear as a function 
of external or internal time delays as has been demonstrated by the Mackey Glass model 
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of blood circulation [18|, the Cheyne Stokes respiration and the forearm tracking with 
visual delayed feedback |HJ. 



The interesting properties of nonlinear delay differential equations have been mainly 
investigated in numerical studies. Therefore it becomes desirable to substantiate these results 
- at least in comparably simple situations - by analytical methods. An interesting result in 
this direction has been recently obtained in 0] by rigorously analyzing the instability of a 
time independent reference state. The application of the theory to a delay induced Hopf 
bifurcation has been confirmed by numerical as well as experimental studies. We note, 
that a different method, which is based upon a multiple scaling analysis, has been recently 
demonstrated in |17[| . 

Here, however, it is our aim to generalize the approach of || by starting from a time 
periodic reference state and by analytically investigating its stability. 



Our paper is organized as follows. In Section II we introduce two methods for approx- 
imately determining a time periodic reference state. Section III then develops its linear 
stability analysis. The resulting Floquet theory leads to a homogeneous vector valued re- 
currence relation determining the Floquet eigenvalues and its corresponding eigensolutions. 
In Section IV we offer two solution methods for this recurrence relation which are based 
on matrix valued continued fractions. Section V completes the Floquet theory by studying 
the adjoint problem. Eventually Section VI is devoted to a short summary and several 
conclusions in view of possible future work. 

For a numerical derivation of the Floquet exponents and the corresponding eigenvectors 
for this case we refer to 



II. DETERMINATION OF THE TIME PERIODIC REFERENCE STATE 

We assume that the dynamical behavior of the system under consideration can be char- 
acterized by a state vector q(t) in an n dimensional state space T and that the underlying 
equation of motion is an autonomous delay differential equation of the general form 

— * 

Here N denotes a nonlinear vector field which depends on the state vector q at the times 
t and t — t, respectively, with r representing the time delay. The set {a{\ describes the 
control parameters which measure external influences on the system. We assume that these 
control parameters are kept fixed so that we can omit them in our notation. 

The treatment of the nonlinear problem ([!]) close to an instability strongly depends on 
the chosen reference state. A theory for an instability of a time independent reference state 
has been recently developed in |J . Here it is our aim to generalize this method to situations 
where we start from a time periodic reference state. As such states cannot be expressed by 
closed analytical forms, it becomes necessary to describe them by using proper approxima- 
tion schemes. In this section we extend methods which have been developed in the realm of 
ordinary differential equations towards delay differential equations. The Poincare Lindstedt 
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approximation allows to determine the time periodic reference state for small values of a 
parameter, whereas its improvement, the Shohat method turns out to possess a wider range 
of applicabilities @,[^| . 



A. The Poincare Lindstedt expansion 



A perturbative approximation method, as the Poincare Lindstedt expansion, relies on 
the existence of a suitable smallness parameter //. Dealing with the nonlinear differential 
equation (|l|), we have to distinguish in general two different origins for such a smallness 
parameter \i. On the one hand, the smallness parameter \i can be generated by a delay 
induced instability. Then it measures the relative deviation of the time delay r from the 
critical value r c above which the delay induced time periodic reference state exists. This 
case occurs, for instance, in the electronic phase locked loop with time delay |J where 
the underlying model equation reveals a Hopf bifurcation at some r c . Considering the 
corresponding normal form [^] 

§ = »Z- 9 |Z|»Z (2) 

we may choose \i = JJr — t c )/t c . On the other hand, the smallness parameter // can also 
coincide with one of the given control parameters of the system. An example is provided 
by a harmonic oscillator with the frequency ujq which is driven by a nonlinear time delayed 
perturbation: 

q"(t,fj) +ulq{t,n) = iif(q(t,n),q'(t,n),q{t-T,n),q'{t-T,y)). (3) 

Here q(t,fi) denotes a scalar variable, the prime abbreviates the derivative with respect to 
the time t and / represents a nonlinear function of its arguments. 



For the sake of simplicity we now discuss the Poincare Lindstedt expansion not for the 
general delay differential equation (P but only for the model equation fl3|). We start with 
the situation of a vanishing smallness parameter // where the solution of (j3|) is a periodic 
reference state 



q(t,0) = q(t + T o ,0), (4) 

with T = 2tt/cu denoting the period of the unperturbed oscillator. Switching on the 
smallness parameter /x, this state will be transformed to a new periodic state, which can be 
described by 

q {t^) = q {t + ^- y ^. (5) 

In the following it becomes useful to explicitly take into account the frequency shift from uj 
to u(fj) by rescaling the time t according to 



3 



£(t) = u{fj)t 



(6) 



Introducing the new variable 

t(£ in — n I 



= 9 ( -T-T./^J > ( 7 ) 



which is 2ir periodic in £, 

x(£,/i) =x(e + 27r,/x), (8) 
we can rewrite the equation of motion @ as: 

w(u) 2 x(£, u) + uj 2 x(£, fi) = fif(x(£, (i),x(£, fi),x(£ - u(ji)t, fi),x(£ - w(ji)t, //)). (9) 
The dot indicates the derivative with respect to the dimensionless new time variable £. 

As already mentioned, we assume that fi represents a small quantity so that we can 
expand the frequency u)(]i) and the periodic orbit x(^,fi) in powers of fi according to 

fi) = x (0 + + u 2 x 2 (0 + • ■ ■ , (10) 
a>(/x) = uj + jjwi + [i 2 LU 2 H . (11) 



In addition to the similar procedure for ordinary differential equations |2^ , [23|] we have also 



to consider a corresponding expansion of the time delayed terms in (|9]). This is achieved by 

x(t ~ w(/i)r, fi) = x (£ - lu t) + fi (xi(£ - lu t) - luitx (£ - lu t)) H , (12) 

~ u{h)t : fi) = x (£ - uj t) + n {x x {i - uj t) - UJlTX (i - UJqT)) H . (13) 

If we apply these expansions to the equation of motion @ and combine terms of the same 
power of /i, we obtain in each order a system of inhomogeneous linear ordinary differential 
equations of second order: 

£o(0+zo(0 = o> ( 14 ) 

£i(0 +^i(0 = -2— x (0 + -j/MOj^oCO^o^ - ^ot),£ (£ - w r)), (15) 
u; c^o 



X n (£)+X„(£)= J n (£). (16) 

The inhomogeneity /„(£) which appears in the nth order (|I^) is purely determined by the 
lower order terms x m (£), < m < n. We have to guarantee that our periodicity condition (H) 
is fulfilled in each order of the perturbation theory. However, if the Fourier expansion of the 
inhomogeneity I n {Q includes multiples of the first harmonic terms which are proportional 
to sin(£) or cos(£), the solution x n (£) of ([TBD contains aperiodic secular terms of the form 
£ sin(£) or £ cos(£), respectively. We can avoid these aperiodic solutions by demanding 

r J n (0 sin(0de = 0, 4(0 cos(0de = 0. (17) 

Jo Jo 
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In order to fulfill these two conditions we need two independent parameters. Here we choose 
the constant uo n as the first parameter, whereas the second one can be chosen by imposing 
suitable initial conditions for x n _i(£)> for example 

x n _ x (0) = A n _ 1; x n _ x (0) = 0. (18) 

In this way we obtain a systematic approximation scheme to determine our time periodic 
reference state order by order for small values of the parameter p. 



B. The Shohat expansion 

In a situation where the parameter p is not a small quantity the Poincare Lindstedt 
expansion for the calculation of the time periodic reference state has to be modified. This 
can be achieved by introducing a new smallness parameter p(p) by the prescription 

M = T~7 (19) 
1 + p 

which maps the interval [0, oo) of p onto the interval [0, 1) of p. The resulting method of 
the Shohat expansion can be described as follows. The equation of motion ([)]) is multiplied 
by p 2 . In doing so we become able to expand the periodic reference state x(£, p) as well as 
the product pui(p) with respect to p and obtain 

p) = X (0 + p(/i)Xx(0 + P (p) 2 X 2 (0 + • • • , (20) 

fjwQi) = p(ii)n + pipfri! + p{p) 3 n 2 + ■■■. (21) 

In order to guarantee that the frequency uj{p) approaches the frequency u) Q of the unper- 
turbed harmonic oscillator in the limit p —>■ we have to choose flo = uj . ^From fl2(f ) and 
(pl|) and the inversion of the relation ( |T9"D 

p = M , (22) 

we deduce the expansions 

uj(p) = n + P (p) (fix - n ) + P (p) 2 (n 2 - no + • • • , (23) 

x(£ - uj(p)r, p) = X (£ - n r) + 

+p(p) [X x (e - n r) - (fix - n )rX (e - ftor)] + • • • • (24) 

The further application of the Shohat method is completely analogous to the Poincare 
Lindstedt approximation scheme. It has been conjectured without proof |22], that the 



method works for arbitrary parameter values p > 0. We thus note, that although to our 
knowledge no known counterexample for this conjecture exists, the vailidity of this expansion 
has to be confirmed for each case individually. 
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III. STABILITY OF THE TIME PERIODIC REFERENCE STATE 



We now generalize the linear stability analysis of a time independent reference state 
developed in to the case of a time periodic reference state. To that end we return to the 
general form of the delay differential equation ([!]) and rescale the time according to 

^q(0 = ±Nm, ftt - wr), {*})■ (25) 

Here uj = uj(fi) abbreviates the frequency of the time periodic reference state, henceforth 
denoted by q°(£) = q°(€,fi). 

Following the original notion of Krasovskii and Hale |5],[7| as well as its detailed elabo- 
ration in || we generalize the n dimensional state space T to an infinite dimensional state 
space C. This allows to embed the given delay differential equation fl25|) in the context of 
functional differential equations. It turns out that this reformulation represents the ade- 
quate framework for a linear stability analysis around a time periodic reference state. The 
resulting Floquet theory leads to a homogeneous vector valued recurrence relation which 
determines the Floquet eigenvalues as well as the corresponding Floquet eigensolutions. 



A. Formulation of the problem in the extended state space 



It appears that solutions of the delay differential equation fl25|) for times £ > depend 
on initial values of the state vector <f(£) in the entire interval [—ujt, 0]. Therefore we have to 
complete (|25|) with the initial condition 



q(8) = g(6), -ujt <6<0, (26) 
where g is a given continuous vector valued function in a suitable function space C. The ini- 



tial value problem (p5|), (26) then maps the function g onto a trajectory in the n dimensional 
state space V. Therefore the problem arises that different initial vector valued functions g 
may yield crossings of the corresponding trajectories in V. This means that the point-wise 
uniqueness of solutions can not be assured when we restrict our considerations to the state 
space T. 



In order to solve this problem one may introduce the extension of the finite dimensional 
state space V to an infinite dimensional function space C where the initial vector valued 
function g is defined. According to Krasovskii and Hale this is achieved by regarding 
the trajectory <f(£) in the original state space V during the time interval [£ — ujt, £] as a single 
point in the extended space C (compare Fig. 1): 

Qk( e ) = + 0), -ujt<9< 0. (27) 

The dynamics of the delay system can then also be described in the extended state space C 
by introducing the nonlinear solution operator 

&(0) = (T(£)0)(0), -ujt<9<0. (28) 
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Its uniqueness is expressed by the fact that the operator T(£) has the properties of a semi 
group, that is 

nz+v)=nzmv), z,v>o, no)=x, (29) 

where I denotes the identity operator. We now have to reformulate the original initial value 
problem Q25D, fl2"B|) in the extended space C. To this end we formally differentiate Q2"8| ) with 
respect to the time £ 

4f e (0) = (A® {9), -ur<9<0. (30) 
Here A denotes the infinitesimal generator which corresponds to the solution operator 

(A® (9) = lim 1 [(T(e)£) (0) - £(0)] • (31) 

By evaluating this limit separately for the interval — lot < 9 < and for the point 9 = we 
obtain the explicit expression 

(Md (0) = l Te m > ~"^ e < °> (32) 
l^(.)], e = o. 

The nonlinear functional M is constructed as follows. We assume that the original vector 
field N in f ^q) can be expanded into powers of its arguments q(£) and <f(£ — cjt). A typical 
term of second order in this expansion has, for instance, the form 

N$ qj (0q k (Z-UT), (33) 

where the explicit components of the respective vectors have been introduced and summation 
is understood over dummy indices. The representation of <f(£) and <f(£ — ujt) can be given 
in terms of the extended state space C by taking into account the relation fl2"7] ) : 

q(0=f° d98(9)q^9), qXZ-UT)=[° d95(9 + ur)q^9). (34) 

J —UJT J — UJT 

If we apply this procedure to every term in the series expansion and collect terms of the 
same order in the extended state vector q^, the nonlinear vector field N becomes a vector 
valued functional M with the components 

M&(0] = E / 9 k )q Ul {9 l )...q Uk {9 k ) 1 (35) 

where the Qu l4 ..j k {9i, . . . , 9 k ) represent matrix valued densities. Thus we have reached our 
first goal, namely to derive a nonlinear functional differential equation for the problem 
formulated in (|25|), (|26|). 
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B. The linearized equation of motion 



According to the prescription Q27D the time periodic reference state q°(£) in the state 
space T transforms into q^(9) in the extended state space C. In order to test its linear 
stability we insert the ansatz 



(36) 



into fl30D , (23). Dropping the tilde we obtain in the linear approximation for the infinitesimal 
deviation q^(9) 



where the linear infinitesimal generator Al becomes explicitly time dependent: 



(37) 



de'n 



-ujt < 9 < 0, 
9 = 0. 



(38) 



The matrix valued density f2§(#) can be written as a functional derivative of M evaluated 
at the time periodic reference state 



n € (0) 



mo) 



(39) 



<?«=<?£ 



C. Transformation of the linear problem 

Due to the fact that the reference state q2(9) is 2% periodic with respect to £, the matrix 
valued density 0^(6') in (^) is time dependent with the same period. We thus perform a 
Fourier expansion of the matrix ft ^(9) : 

oo 

n f (0)= £ n k {9)e ik *. (40) 

k=— oo 

In close analogy to the Floquet theorem for ordinary differential equations p5| we try to 
solve (|3)-(|0D by the ansatz 

m = e*%(9). (41) 

Here A denotes the Floquet eigenvalue and <fif(9) = <^ +27 r(#) is a 2vr periodic Floquet eigen- 
solution for which we also perform a Fourier expansion: 



(0) = e in{oy ntL - (42) 
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— # 

In order to determine the Fourier components A (#) we insert the hypothesis (PJ), (H) into 
the linearized equation of motion (j37p-(|40D. We now have to consider separately the interval 
—lot < 9 < and the point 9 = 0. In the interval —lot < 9 < we conclude that the Fourier 
component A (#) has the form 

fc(0) = ^J x+ ^ e . (43) 

For the case 9 = we find 

OO OO OO „Q 

E A (A + m)e (A+iri) « = EE/ ^^fc(^)e (A+in)e e (i(fc+n)+A) «^. (44) 

n=— oo n=— oo k=— oo — tJr 

We now introduce the matrix valued quantity 



J k,n 



° dflft fc (#)e (A+in)e (45) 



and a new index ra(ra) = n + k. Comparing the contributions of the various Fourier com- 
ponents in (H4T) and dropping the tilde, we obtain a homogeneous vector valued recurrence 
relation for the Fourier components A : 

oo 

0= E [U,n-k-5 k ,o(X + tn)I\$ x n _ k . (46) 



Thus we are left with the problem to construct an approximate solution to ([46]) which leads 
to both the Floquet eigenvalues A and the corresponding Floquet eigensolutions. 



D. Remark 



In the Floquet theory of ordinary differential equations |25| it is shown that the derivative 
of the time periodic reference state represents a Floquet eigensolution where the real part 
of the corresponding Floquet eigenvalue vanishes. This statement remains valid for delay 
differential equations as can be seen as follows. As the time periodic reference state q^(6) 
satisfies the nonlinear equation of motion (|30|), (p^), a differentiation with respect to the 
time £ leads to 



d dqf(6) 



d dqj?(6) 



d9 



di 

d0' 



-LOT < 6 < 0, 



9 = 0. 



(47) 



A comparison with d37|) (p9|) reveals that the derivative of the time periodic reference state 
qf(9) indeed fulfills the linear problem. Due to ( f4~T| ) and (|42| ) it therefore possesses the 
general form 



d£ 



E 



(48) 



^From the 2tt periodicity of q®{9) and its derivative 
of its Floquet eigenvalue A has to vanish. 



§D we the conclude that the real part 
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IV. MATRIX VALUED CONTINUED FRACTIONS 



We consider two methods which enable us to approximately solve (EBD for the Floquet 
eigenvalues A and for the Fourier components 0^ of the Floquet eigensolutions. In the first 
part we formulate a new method based on n diagonal continued fractions. In the second 
part we show that this solution method is equivalent to a formulation with tridiagonal 
continued fractions introduced by Risken PB| . It turns out, however, that the first method 
is much simpler to handle as the necessary inversion of matrices can be performed in a 
low dimensional space. Furthermore the criterion for truncating higher order terms in the 
smallness parameter \i can be formulated more precisely in the framework of the first method. 



A. Pentadiagonal recurrence relations 



We apply the method of matrix valued continued fractions in order to solve the vector 
valued recurrence relation ([16]) approximately. In order to avoid overloading the notation 
we restrict ourselves for the time being to the pentadiagonal case where the summation in 
(ISP is performed for — 2 < k < 2: 

= L_ 2 ,n+20n+2 + L -l,n+10n+l + [ L 0,n ~ (A + «n)I]<^ + Li^-i^.j + L 2 , n -2^n-2- (49) 

We start with defining a set of ladder operators S™ for m = ±1, ±2 which relate neighboring 
Fourier components via 



n+m 



(50) 



This definition implies the following useful relations between different ladder operators: 

-i 



c+i 

°n-l 



o+l O+l 

°n.+l°n. 



2+2 



(51) 



Applying the definition (|50f) of the ladder operators, the pentadiagonal recurrence relation 
(03T) can be rewritten as: 







L_ 2in+2 S+" + L^n+iSiJ; 1 + [L , n - (A + m)I] + L ljn _xS 



-i 



(52) 



We now express the ladder operators S™(m = ±1,±2) in terms of the matrices Lfc in as 
well as the operators S™ ±1 , S™ ±2 . In order to evaluate this dependence, we isolate the term 
). Then the equation assumes the form 



n+lVn+1 



[Lo,n — (A + in) I] + L_ 2in+2 S^ 2 + L lin _ 1 S n 1 + L 2in _ 2 S n ' 



(53) 



Shifting the index from nton-1 and applying the definition (fi^-i 
the validity for all ^ the operator relation 



2-1ZX 



we obtain from 



J 0,n- 



-i-(A + i(n-l))I| 



-2,n 



+iS^! + L 1]n _ 2 S n i 1 + L 2in _ 3 S„ 2 1 L_ 1)n . (54) 



Similarly we construct the operator relations 
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Q+l 

s- 2 



[Lo,n+i - (A + i(n + 1))I] + L_ 2 ,„+ 3 S^f j + L_ 1>n+2 S+^ + La.n-iS^ L 1>n , (55) 



[L„,„-2 - (A + i(n - 2))I] + L-Ln-iS+Ia + L lin _ 3 S 



3-1 
n-2 



+ 



3-2 



[Lo,n+ 2 - (A + i(n + 2))I] + L_ 1)n+3 S+| 2 + L li7l+1 S-j 



+2 



T 



We perform an iteration procedure by starting from 



for the case n — 0, 







L-2,2S([ 2 + L-ijS^ 1 + [Lq,o — AI] + Li 5 _iS 1 + L^^Sg ' 



L-2,n, (56) 

_1 L 2 ,„. (57) 



(58) 



and by recursively inserting the recurrence relations of the ladder operators (p^)-(|57j) . Writ- 
ing the successive inversions of the matrices formally as fractions, we may visualize this 
iteration procedure by a schematic representation of a pentadiagonal matrix valued contin- 
ued fraction. This is illustrated in Fig. 2 where each term is represented by a horizontal line 
and where bifurcations are related to the terms of the ladder operators S™(m = ±1, ±2). 



This far we discussed the solution of the vector valued recurrence relation ([52]) in the 
pentadiagonal case for m — ±1, ±2. However, our method can be correspondingly extended 
to the general case where all Fourier components <^ are coupled to each other. To this 
end we introduce ladder operators S™ with arbitrary m according to (0), where we identify 
S° = I. The homogeneous vector valued recurrence relation ( [46"D then yields a corresponding 
one for the ladder operators S™: 

oo 

0= [U,n-k-h fi (\ + tn)l]S- k . (59) 

k=— oo 

An iteration procedure similar to (|51])-(|5"7j) finally leads to a homogeneous equation for the 
Fourier component 0q, 

M(A)$ = 0, (60) 

where the resulting matrix M(A) consists of an infinite number of matrix valued continued 
fractions transcendentally depending on the Floquet eigenvalues A. Therefore the Floquet 
eigenvalues A are determined from the condition that the determinant of the matrix M(A) 
vanishes: 



detM(A) = 0. (61) 

Once the Floquet eigenvalues A are known, the yet unknown Fourier component 0q * s deter- 
mined up to a constant from solving the homogeneous equation (j60|) . All Fourier components 
<fin of the Floquet eigensolutions are then calculated by successively applying the ladder op- 
erators S™ starting with (f)^: 

& = S£4 (62) 

Note that the remaining normalization constant in 0q has to be fixed by an adequate 
biorthonormality condition which will be discussed in Section V.D. 
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In applications, however, it is impossible to exactly evaluate the infinite number of matrix 
valued continued fractions. ^From an analytical point of view we can therefore expect that 
this solution method allows at most to approximately determine Floquet eigenvalues and 
the corresponding eigensolutions. To this end we recall that the starting point of our linear 
stability analysis, i.e. the time periodic reference state, is only known as a finite power series 
in the smallness parameter /i. As a consequence the whole calculation can be simplified by 
approximately neglecting higher order terms in the smallness parameter /x. In particular it 
becomes sufficient to restrict the vector valued recurrence relation fl46|) to a finite number of 
terms, so that the subsequent iteration procedure only leads to finite number of matrix valued 
continued fractions. Furthermore each continued fraction can be evaluated in the leading 
order of the smallness parameter fi. In spite of these successive expansions, the continued 
fractions have the property that the approximate results rapidly converge towards the exact 
values if the leading order in the smallness parameter /j is increased. 



B. Sketch of Risken's tridiagonal formulation 



Following Risken |26[ we show that the n diagonal matrix valued continued fractions can 
always be cast into a tridiagonal form. For the sake of simplicity we demonstrate this only 
for the pentadiagonal recurrence relation (f49|), but the general case is treated along similar 
lines. We start with distinguishing between even and odd indices n in the pentadiagonal 
recurrence relation (E9|): 



— L_ 2 ,2n+202n+2 + ^-l,2n+X<t>2n+l + 



+ [L ,2„ - (A + l2n)l}4> X 2n + + L2,2„-202n- 2) (63) 



— L-2,2n+302n+3 + L-l,2n+202n.+2 + 



+ [L ,2n+l " (A + l{2n + l))I]02n+l + Ll,2n02n + U,2n-A-V (64) 

We now construct new vectors according to the prescription 






<i = ST > *» = rT , K-i = Vr • (65) 



Additionally we define the matrices 

Q__ I L_ 2 ,2n+20 | r» ( I J 2 ) 2n-2l J l I 2n-l \ 

-l,n+l — x t 5 ^l,n-l nT ; 

\ -Li-l,2n+2-^-2,2n+3 / \ Ulj2,2n-1 / 

n ( [L 0) 2n - (A + 2m)I]L_i, 2 „ + l \ 

^ n { L 1)2n [L 0) 2n+i - (A + i{2n + 1))I] ) ' {m > 

Due to these definitions both pentadiagonal recurrence relations (|63D, (|6^) can be combined 
in the following way: 

Q_ lin+1 ^ +1 + Q , n ^ + Q +1 , n _ 1 ^_ 1 = 0. (67) 

Thus we have reached our goal to find a tridiagonal vector valued recurrence relation. At 
this stage we define again ladder operators with the property 
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$ A „ = R+$ A 



n-1 



R~$ A 



(68) 



These ladder operators can be determined when we rewrite the tridiagonal recurrence rela- 
tion (|67|) as 







Q-l,n+lR+ + Q ,n $n + Ql,«-1^-1 



and, similarly, as 



— Q-l.n+l^n+i + Qo,n + Ql,n-lR,,, 



(69) 



(70) 



Comparing these results with the original definitions fl6"8|) and shifting the index we find 
relations for the ladder operators themselves: 



Q±l,n=F2R-nTl + Qo,n=Fl Q=Fl,n- 



(71) 



Again it is sufficient for our purpose to solve the tridiagonal recurrence relation (|67j) for the 

case n = 0: 







Q-i,iR-o" + Qo,o + Qij-iR-o 



(72) 



Repeated application of the operator relations ( [HP then yields a tridiagonal matrix valued 
continued fraction. A schematic representation of the iteration method is presented in Fig. 
3. 



V. FORMULATION OF THE ADJOINT PROBLEM 

In general, the linear infinitesimal generator Al is not self adjoint in the extended state 
space C. Therefore it becomes necessary to define another extended state space & dual to 
C and to investigate the properties of the adjoint infinitesimal generator A L . In order to 
relate the linearized problem with its adjoint it turns out that the canonical bilinear form for 
ordinary differential equations is not appropriate. In the case of delay differential equations 
a modified bilinear form has to be introduced. 



A. The bilinear form 



The choice of the canonical bilinear form for delay differential equations is motivated 



from the Fredholm alternative. To this end we consider the inhomogeneous version of ( 3~7|) 



Al - 



d 



ok (0) = -ut<9<o 



(73) 



with a 27r periodic vector valued function x^{9) = Xt+2n(@)- We try to construct a particular 
solution q^{6) of ( [73] ) by the Floquet ansatz flflD with 4>^{6) = 0^ +2ff (6')- Inserting the Fourier 
expansion ( f42| ) for <fi^{0) and a corresponding one for the inhomogeneity x^{9), we obtain 
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E ([A L -\- in) fa) E UO)e mi . (74) 

n=— oo n=— oo 

Taking into account the definition (|38| ) of the infinitesimal generator Al in the interval 
— lot < 9 < 0, we conclude from ([74]) the general form of the Fourier component 0^(#): 



= <^(0)e( A+m ) e + / dse^ m ^Us)- (75) 



Correspondingly ( |T4] ) determines for the point 8 = the yet unknown initial condition 0^(0). 
With the definition (flop it fulfills an inhomogeneous vector valued recurrence relation: 



E t L M-fe - 4,o(A + m)I] <_ fe (0) 

k=— oo 



Xn(0) - E / W ^e( A+J ("- fe ))^- s ^ fc (^)xn- fc (s). (76) 

fc=-oo J -^ T J ° 



This result suggests how to introduce both the dual space and the bilinear form. We 
assume that consists of n dimensional vector valued functions defined on the interval 
[0, lot] and that the bilinear form is given by 

U e ))^ = (^(0), <M0)) - d6 J° ds - 9), Q^s-emd*)) (77) 

for all (f)£ G C and $ G C^, where <, > denotes the usual canonical scalar product. Note 
that each delay system and each time periodic reference state induces its own bilinear form 
due to (|3^). Furthermore we observe that the explicit time dependent bilinear form (|77|) 
for a time periodic reference state reduces to the corresponding one for a time independent 
reference state ||. 

With the bilinear form (|7"T| ) the inhomogeneous recurrence relation ( |T6| ) can be rewritten 
according to 

OO 1 ^27T 

E [U,n-k - 4,o(A + m)I] <gU(0) = — J o d£ (A[ n (s), %(ej) 6 , (78) 



fc=— oo 



where the matrix valued functions A^ n (s) are given by 

A\ in (s) = e-^e-^+^'I, 0<s<ut. (79) 

Thus we obtain the following Fredholm alternative for solving the inhomogeneous equation 
(f73|). If the parameter A does not coincide with a Floquet eigenvalue, we read off from 
([78D that there exists a unique solution. Otherwise we can only expect a solution if the 
inhomogeneity \£ fulfills a solvability condition which involves the bilinear form (|77|) . This 
solvability condition will be concretized below after having defined the adjoint operator A L 
and its corresponding Floquet eigensolutions. 
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B. The adjoint operator 

The bilinear form ( [77p can be applied to describe the evolution of the linearized delay 
system also in the dual extended state space C'. To this end we require that the bilinear 
form between the state vector q^ 6 C and its dual q*g G becomes time independent: 

= ^(ql(s),m)^ (80) 

As the bilinear form (|77| ) does explicitly depend on the time £ via the matrix valued density 
0,^(6), we derive from (|37|) and (^) the evolution equation in 

^(s) = -(Atqp)(s), 0<s<ujt, (81) 

where the adjoint infinitesimal generator A L obeys: 

(Alqlq^ = (qlA L q^- J°^d6 jf* ds (q^s - 9), ^+^(9)^) V (82) 

When we use the definition ([38]) of the infinitesimal operator Al we obtain from (|82|) after 
a partial integration the following expression for the adjoint infinitesimal generator A\: 



d 



s), < S < UJT, 



'y ds'ql(s')^ +s/ (-s'), 8 = 0. 



C. The adjoint recurrence relation 

We are now in the position to solve the adjoint problem defined by (^l|) and (0). In 
close analogy to the procedure in Section HI.C we perform the Floquet ansatz 

fl?(«) = e- x ^l\s) (84) 

with the 2tt periodic adjoint Floquet eigensolution 

$\s) = Y,$\s)e-^. (85) 

3 

Evaluating (jSl~l) and (^) in the interval < s < ujt fixes the form of the Fourier components 
according to 

$\s) = $V< A+i '> (86) 

whereas the case s = leads to the corresponding homogeneous vector valued recurrence 
relation 
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0= E [L fcj - MA + *J)I] • (87) 

k=— oo 

In order to solve fl8~7|) for the Fourier components tp] x of the adjoint Floquet eigensolutions 



3 

and the respective Floquet eigenvalues A, we proceed along similar lines as in Section IV. A. 
First we define adjoint ladder operators Z™ with arbitrary m and the identity Z° = I in 
analogy to (|50|) : 

$ x +m = $ X Z™ . (88) 



Inserting (88) in the homogeneous vector valued recurrence relation (87) we then obtain a 



corresponding one for the adjoint ladder operators: 

oo 

0= £ Z*[L k}j -5 k)0 (\ + ij)l}. (89) 

k=— oo 

A careful comparison between (^) and (|89|) reveals that the recurrence relations for the 
ladder operators S™ and their adjoint Z™ are not independent from each other. Indeed, 
they are mapped into each other by the prescription 

Lfc,n-fcS n fc = Z n fc L_fc in . (90) 

This means that the adjoint ladder operators Z™ can be immediately calculated, once the 
ladder operators S™ are known. However, this does not imply that the solution of the 
adjoint problem directly follows from the linear problem. Iteratively inserting the operator 
recurrence relation (|39"D in the vector valued recurrence relation (|87|) for j = yields with 

m 

4 X M(X) = 0. (91) 



Thus the adjoint problem leads to the same condition (|61[ ) for the Floquet eigenvalues A, but 
the Fourier component ipQ X has to be determined independently from the Fourier component 



% defined by © 




With these definitions we are now able to concretize the Fredholm condition for solving 
the inhomogeneous equation (|73|). Multiplying ( |78|) from the left with ip^ x , performing the 
summation over all n and taking into account the homogeneous vector valued recurrence 
relation (0), we yield 

E ^ X AUs),M0)) =0- (92) 

-OO / £ 

Due to ([79]), (|H3|), ( |8"6j ) this solvability condition takes the concise form 

hr*$ x ' m \ = °- (93) 

Only in the special case that all quantities do not explicitly depend on the time £ this 
reduces to the usual Fredholm condition, i.e. the inhomogeneity \(, must be orthogonal to 
the respective eigensolution of the adjoint operator A L f|. 
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D. Biorthonormality relations 

Finally we show that the Floquet eigensolutions 0| and ip^ x of the infinitesimal generator 
Al and its dual A L , respectively, can be chosen to form a biorthonormal set in C and C^. 
First we derive the biorthogonality condition for different Floquet eigenvalues/i 7^ A. Using 
the explicit expression for the bilinear form ([77]) and applying our previous results PDp- 



and (|85f), (|86|) we obtain 



00 



00 „q „g 

- J2 e i{n+k -^ J d9 I (^) A , ^ fc (6i)^) e (^(i-fc))« e (/,-A+i(^+fc-i)) s ^_ (94) 

k=~ 00 







An integration with respect to s yields for the second term on the right hand side 

00 o p {n+in)6 _ J\+i(j-k))9 

- ' L de w ,-A + ,- (n+t - jr < 95 > 

With the definition (flop of the matrices L fc n this reduces to 



.f" u-\ + i(n + k-j) 

n,3,k=— 00 " \ J / 

so that we yield from the homogeneous vector valued recurrence relations (|46|) , fl87|) 



^From (p4[) (p7|) we conclude the biorthogonality for /1 7^ A 

>« 



- £ (97) 



tA ^),^W) =0. (98) 



5 v^rg 

To normalize the biorthogonal set of Floquet eigenf unctions, we introduce a proper normal- 
ization constant in a symmetric way: 

& = N x $l ^ = N x & n \ (99) 

^From the requirement 

($ A ( S ),$(0)) e = l (100) 

we then determine the normalization constant N\ by performing similar calculations as 
above: 

1/2 



(101) 

\ \ ./ 1, IT / / 

n ,]=— 00 

As expected the normalization constant iV\ does not explicitly depend on the time £. Sum- 
marizing the results fl9"8|) and ( |100| ), the biorthonormality relation reads 



c ^),#(*)), = <W- (102) 
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VI. SUMMARY AND CONCLUSIONS 



The present paper was devoted to systematically develop a Floquet theory for delay 
differential equations. At first we approximately determined a time periodic reference state 
by extending two standard methods for ordinary differential equations, namely the Poincare 
Lindstedt and the Shohat expansion. Then we tested the stability of this reference state 
by constructing Floquet eigensolutions and their corresponding eigenvalues from matrix val- 
ued continued fractions. Finally the Floquet theory was completed by studying the adjoint 
problem. The applicability of our Floquet theory was demonstrated in In particular 



our analytical treatment provides means of understanding the mechanism of the continuous 
control of chaos by self controlling feedback |[7|,[2]|. Previous investigations have indicated 



that it becomes crucial to decide whether an observed stabilized limit cycle corresponds to 
an unstable cycle of the system or is produced by the control mechanism itself p9|.|30|. 



As the Floquet theory represents a linear stability analysis for a time periodic reference 
state there still remains the nonlinear problem to construct the normal form for an emerg- 
ing instability. We expect that this problem can be tackled in a similar way as in |J where 
synergetic methods p4| , p5| are extended to investigate delay differential equations in the 



local neighborhood of a time independent reference state. Also close to the instability of a 
time periodic reference state the inherent time scale hierarchy should allow to adiabatically 
eliminate the fast modes by using projectors which are induced by the bilinear form ([77|) 
of the linear stability analysis. As in || the resulting order parameter equations for the 
slow modes should turn out to be of the form of ordinary differential equations. We stress 
that the normal form theory is indispensable for classifying the instabilities of time periodic 
reference states. Whereas the linear stability analysis is sufficient to identify the instabilities 
of time independent reference states, this is no longer true for time periodic ones 
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